Tentative evidence for excess GeV-scale gamma rays from the galactic center has been corroborated by several groups, including the Fermi collaboration, on whose data the observation is based. Dark matter annihilation into standard model particles has been shown to give a good fit to the signal for a variety of final state particles, but generic models are inconsistent with constraints from direct detection. Models where the dark matter annihilates to mediators that subsequently decay are less constrained. We perform global fits of such models to recent data, allowing branching fractions to all possible fermionic final states to vary. The best fit models, including constraints from the AMS-02 experiment (and also antiproton ratio), require branching primarily to muons, with a $\sim 10-20\%$ admixture of $b$ quarks, and no other species. This suggests models in which there are two scalar mediators that mix with the Higgs, and have masses consistent with such a decay pattern. The scalar that decays to $\mu^+\mu^-$ must therefore be lighter than $2m_\tau \cong 3.6$ GeV. Such a small mass can cause Sommerfeld enhancement, which is useful to explain why the best-fit annihilation cross section is larger than the value needed for a thermal relic density. For light mediator masses $\sim (0.2-2)$ GeV, it can also naturally lead to elastic DM self-interactions at the right level for addressing discrepancies in small structure formation as predicted by collisionless cold dark matter.
18 pages, 14 figures; v2: updated CMB constraint and added references; v3 corrected direct detection cross section